In this document I compare gene-snp pairs with strong heterogeneity:

knitr::opts_chunk$set(cache=TRUE)
covmat=readRDS("../withoutbrain/covmatnobrain.rds")
pis=readRDS("../withoutbrain/pisnobrain.rds")$pihat
z.stat=read.table("../Data/maxz.txt")
b.stat=read.table("../Data/maxbetahats.txt")
pi.mat=matrix(pis,ncol=44,byrow=T)
lfsr.mash=read.table("../withoutbrain/nobrainlfsr.txt")[,-1]
pm.mash=read.table("../withoutbrain//nobrainposterior.means.txt")[,-1]
pm.mash.brain=read.table("../Data/Aug13withEDposterior.means.txt")[,-1]

standard.error=read.table("../Data/standard.error.txt")
pm.mash.beta.brain=pm.mash.brain*standard.error

pm.mash.beta=pm.mash*standard.error[,-c(7:16)]

pm.ash=read.table("../Comparison with Univariate Ash/univariate.ash.pm.txt")[,-c(7:16)]
lfsr.ash=read.table("../Comparison with Univariate Ash/univariate.ash.lfsr.txt")[,-c(7:16)]
pm.ash.beta=pm.ash*standard.error[,-c(7:16)]

colnames(lfsr.mash)=colnames(pm.mash)=colnames(pm.mash.beta)=colnames(z.stat)[-c(7:16)]


het.norm=function(effectsize){
  t(apply(effectsize,1,function(x){
  x/x[which.max(abs(x))]
}))}


het.func=function(normdat,threshold){
  apply(abs(normdat),1,function(x){sum(x>threshold)})}



matrix.ash.norm=het.norm(effectsize = pm.mash.beta)
uni.ash.norm=het.norm(effectsize = pm.ash.beta)
matrix.ash.norm.brain=het.norm(effectsize = pm.mash.beta.brain)

het.genes.no.brain=het.func(matrix.ash.norm,0.5)
het.genes.with.brain=het.func(matrix.ash.norm.brain,0.5)

Let’s looks at how heterogeneity behaves:

library(qtlcharts)
iplotCurves(curveMatrix = pm.mash.brain[het.genes.with.brain>40,])
## Set screen size to height=700 x width=1000

iplotCurves(curveMatrix = z.stat[het.genes.with.brain>40,])

Now, we look at singletons:

iplotCurves(curveMatrix = pm.mash.brain[het.genes.with.brain<2,],chartOpts=list(curves_xlab="Tissue", curves_ylab="PMMASHBrain"))

iplotCurves(curveMatrix = z.stat[het.genes.with.brain<2,],chartOpts=list(curves_xlab="Tissue", curves_ylab="Z.statBrain"))

Let’s consider what happens when we remove brains:

iplotCurves(curveMatrix = pm.mash.beta[het.genes.no.brain >40,])

iplotCurves(curveMatrix = b.stat[het.genes.no.brain>40,-c(7:16)])

Now, we look at singletons:

iplotCurves(curveMatrix = pm.mash.beta[het.genes.no.brain <2,])

iplotCurves(curveMatrix = b.stat[het.genes.no.brain<2,-c(7:16)])